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Abstract 

In recent years, intensive effort has gone into developing numerical tools for exact quantum me- 
chanical calculations that are based on Bohmian mechanics. As part of this effort we have recently 
developed as alternative formulation of Bohmian mechanics in which the quantum action, S, is 
taken to be complex [JCP 125, 231103 (2006)]. In the alternative formulation there is a significant 
reduction in the magnitude of the quantum force as compared with the conventional Bohmian 
formulation, at the price of propagating complex trajectories. In this paper we show that Bohmian 
mechanics with complex action is able to overcome the main computational limitation of conven- 
tional Bohmian methods — the propagation of wavefunctions once nodes set in. In the vicinity 
of nodes, the quantum force in conventional Bohmian formulations exhibits rapid oscillations that 
pose severe difficulties for existing numerical schemes. We show that within complex Bohmian 
mechanics, multiple complex initial conditions can lead to the same real final position, allowing 
for the description of nodes as a sum of the contribution from two or more crossing trajectories. 
The idea is illustrated on the reflection amplitude from a one-dimensional Eckart barrier. We 
believe that trajectory crossing, although in contradiction to the conventional Bohmian trajectory 
interpretation, provides an important new tool for dealing with the nodal problem in Bohmian 
methods. 

PACS numbers: 
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The challenge of performing quantum mechanical calculations on systems of many degrees 
of freedom has been a major focus of the theoretical chemistry community for almost four 
decades. Since classical mechanics can be applied to systems with tens of thousands of 
degrees of freedoms it is only natural that much attention has focused on methods capable 
of describing quantum effects but requiring only the propagation of classical trajectories. 
One approach that has shown significant p ro g ress in recent years is the use of Bohmian 



mechanics 



(BM)[ii,y,y,Q,H,y,0,y,y, 
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12]. Originally developed in the 1950's as 



a causal formulation of quantum mechanics (QM), in this formulation trajectories evolve in 
the presence of the usual Newtonian force plus an additional quantum force [13|]- Although 
the quantum force is nonlocal, the formulation of BM in terms of trajectories suggests the 
possibility of achieving an efficiency that is compatible with classical trajectory methods. 
Despite its successes, BM suffers from several drawbacks that have prevented it from becom- 
ing an effective numerical tool to date, the most significant of these being the nodal problem 
- a numerical instability in regions where the wavefunction oscillates that eventually leads 
to a breakdown of the numerical scheme 121 ] . 

In this paper we address the nodal problem using a recently developed variation on BM 
that we call Bohmian mechanics with complex action (BOMCA) Q, 13]' where the phase 
S, is taken to be complex. Like conventional BM, BOMCA is formally identical to the exact 
Schrodinger equation; however the quantum force is significantly smaller and more localized 
than in conventional Bohmian mechanics. This comes at the expense of the trajectories 
being complex. 

In our previous publication 14] we focused on the transmitted part of a wavepacket scat- 
tered from an Eckart barrier; we demonstrated tunneling probabilities that were in virtually 
perfect agreement with the exact quantum mechanics down to 10 -7 . Here we complete the 
picture by calculating the reflected part of the wavefunction. Whereas the transmitted part 
is nodeless, the reflected part, for a wide range of energies, has an oscillatory structure. We 
show that the oscillatory structure is obtained automatically in BOMCA simply by including 
the contributions from multiple initial conditions that lead to the same final position. 

The origin of the nodal problem in BM can be traced back to the hydrodynamic equations 
of QM, which constitute the first step in the derivation of the Bohmian formulation. In 
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conventional BM these equations take the form: 

A, + -AA + ^-AS„ = 0, (2) 
m 2m 

where V(x, t) is the potential of the system, m is the mass of the particle, h is Planck's 
constant divided by 2n and the subscripts denote partial derivatives. 

is referred as the "quantum potential". As can be seen from the expression for Q, the 
quantum potential diverges at nodal regions of the wavefunction. Numerically the difficulty 
is even more severe — well before a node is formed, when the amplitude of the wavefunction 
exhibits only nodeless ripples, the quantum trajectories are highly unstable due to rapid 
oscillations in the quantum potential 1^] . 

Since nodes in quantum mechanics arise from interfering amplitudes, it is only natural 
to attempt to solve the nodal problem by applying the superposition principle — to de- 
compose the wavefunction into two nodeless parts and to propagate each part separately 
using trajectories. Indeed, two such methods have been developed, the Counter Propagating 
Wave Method (CPWM)0| and the Covering Function Method (CFM)Q- However, since 
nodeless wavefunctions do not stay nodeless for a general potential, these methods gener- 
ally require a series of time-dependent decompositions of the total wavefunction, which is 
numerically inconvenient and largely arbitrary. 

A somewhat more natural strategy to solve the nodal problem is to apply the superposi- 
tion principle directly to the contribution of the quantum trajectories. The superposition of 
contributions from multiple trajectories is a central concept in the semiclassical literature [it]]; 
however, in conventional BM, the crossing of trajectories in configuration space is strictly 
prohibited. Indeed, the "no-crossing" rule plays a central role in conventional BM: if tra- 
jectories could cross it would undermine the Bohmian interpretation of QM, in which the 
quantum trajectories are candidates for the actual trajectory on which a particle propagates 
in space. Since the BOMCA formulation yields generally an approximation to QM, trajec- 
tories are allowed to cross. In this paper we demonstrate how combining the contributions 
of crossing trajectories yields accurate interference patters. 
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The starting point of the BOMCA formulation [ljjl is the insertion of the ansatz 17l.ll8l.ll9. 



201 ] ip(x,t) = exp [tS'Oe, t)\ in the time-dependent Schrodinger equation, where we allow the 



phase to be complex. This yields a single quantum complex Hamilton- Jacobi equation 
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In the spirit of conventional BM our aim is to solve eq. (T4|) in the Lagrangian approach, that 
is along quantum trajectories. A quantum trajectory is defined by 

dx 1 

— = v(x,t); v(x,t) = — S x (x,t). (5) 

Due to the definition of x as time-dependent in eq.([5]), we write the solutions of this equation 
as x(t;x ) where Xq is the starting point of the trajectory. Unlike conventional BM, the 
complex value of S yields quantum trajectories x(t; xq) that evolve in the complex plane. A 
Newtonian-like equation of motion for v(x,t) is obtained by taking a spatial derivative of 
eq.(H]) and applying eq.([5]); after a short manipulation we obtain 

dv[x(t;x ),t] V x ih 

; = 1 v xx , (6) 

at 

F c /m F q /m 

where we identify F c , F q as the classical and the quantum force respectively. The presence 
of v xx in the quantum force term prevents the first equation in ([5]) and eq.(E]) from being a 
closed set. 

We tackle the problem of calculating v xx along a trajectory by taking iterated spatial 
partial derivatives of eq.([6]). After a short manipulation the result can be written as 

dv {n) y(n+l) %h 

—— = + —v {n+2) -g n , n = 0,...oo, (7) 

at m 2m 

where go = and g n = 2~2j=i i^v^v^^ 1 ' for n > 1. Here denotes the n th spatial 
derivative of v. A similar procedure was used in Refs. in conventional BM. The 

set of eqs.([71) and the first equation in ([5]) are now an infinite but closed set that describes 
a fully local complex quantum trajectory. We may obtain a numerical approximation by 
truncating the infinite set at some n = N, thus replacing eq.Q with a system of iV + 1 
coupled ODEs. Since each equation of motion for v ^ in ([7]) depends on the subsequent 
v (n+2) ^ truncation is done by setting v^ N+1 ^ = v^ N+2 ^ = 0. The initial conditions for the 
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UW'g are gj ven |)y 

^ (n) (0;x ) 



1 d n S x (x,0) 



X = XQ 



—ih 



^(x,0) 



(8) 



m <9x n 

where we have applied the definition from (jSJ) together with 5* = — ihln if), the latter following 
from the original ansatz. xq is the initial position of an arbitrary single trajectory, which is 
generally complex. The equation of motion for the action along a trajectory is similar to its 
classical counterpart, with the addition of the quantum potential: 



dS[x(t; Xo),t] 
dt 



S t + vS x 



—mv 
2 



ih 

V+-v x . 



(9) 



Solving for S and inserting the result into the original ansatz yields the wavefunction 
i/j[x(t] x ),t] = exp {j:S[x(t; x ),t]} at position x(t; x ) in the complex plane. In Ref. [14]], we 
discussed the practical issue of finding trajectories that end up on the real axis x(tf] Xq) G M 
at tf. Once these trajectories are found, the wavefunction on the real axis if)[x(tf]Xo),tf\ is 
reconstructed. 

As a numerical example, consider the one-dimensional scattering of an initial Gaus- 
sian wavepacket ip(x,0) = (2a/rc) 1/ ' i e\-'~ a( - x ~ Xc ^ + ^ Pc( - x ~ Xc ^ from an Eckart potential V(x) = 
D/ cosh 2 (/5x). We take x c = -0.7, a = 30tt, D = 40, /3 = 4.32 and m = 30 (all units 
are atomic units). The system is the same as in our previous publication 14|. The average 
translational energy of the initial Gaussian is taken as E = p 2 c /m = 10 < D. We focus on 
trajectories that end up at a final time tf — 0.995 with x(tf) < 0, and thus contribute to 
the reflected part of the wavefunction. tf is chosen as sufficiently long for the wavepacket 
to scatter from the barrier and interference effects to appear. First we focus on the N = 1 
BOMCA approximation, for which the equations of motion are: 

dx dv V r dv r V*, 

— — v. 
dt 



dt 



' xx 2 



m 



(10) 



with the auxiliary eq. 



dt m 

1 trajectories obey Newton's equations of 



Note that the N 

motion. As indicated in Ref. jl4 | , eqs. fllQI) are closely related to Generalized Gaussian 
Wavepacket Dynamics (GGWPD)[21!|, which also uses complex trajectories. As a result, 
many of the insights in Ref. 2lJ concerning multiple root trajectories can be carried over to 
the case of complex Bohmian mechanics, although we have found that the structure and the 
number of the root branches generally depend on the value of N, the order of truncation. 

In figjl]we plot three branches that contribute to the reflected wavefunction. A branch 
is defined as the locus of initial positions of trajectories that end at time tf at real Xf with 
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FIG. 1: Scattering of a Gaussian from an Eckart potential barrier within the N=l BOMCA 
approximation (the parameters are given in the text). Three branches and six sample complex 
trajectories are depicted. The branches are the locus of initial positions of trajectories that end at 
time tt = 0.995 at real Xf, —1 < Xf < —0.05, and thus contribute to the reflected wavefunction. 
The figure shows two sample trajectories emerging from each branch, one that ends at Xf = — 1 
and one that ends at Xf = —0.05. 

Xf corresponding to a reflected segment of the wavefunction, Xf G [—1, —0.05]. Two sample 
trajectories are depicted emerging from each branch, one that ends at position Xf = — 1 
and one that ends at Xf — —0.05. Thus, to each final position (for this segment of the 
reflected wavefunction) there correspond three initial positions — one originating from each 
of the three branches — and the wavefunction should therefore include contributions from 
all three branches. As in Ref. 21] ; at short times only one branch contributes to the final 
wavefunction. We refer to this branch as the "real branch" since it incorporates a trajectory 
that stays on the real axis at all times (the trajectory that initiates from xq = x c ). At longer 
times, secondary branches begin to make a significate contribution to the final wavefunction. 
There is apparently no fundamental limitation on the number of the secondary branches, 
although in the section of the complex plane depicted in figJUno other branches were found. 
In ng{TJ branch (2) is the real branch while branches (1) and (3) are secondary branches. 

In figJ2ta) we plot \ipj(x,tf)\ = \ exp[iSj(x, tf)/h)\, j = 1,2,3, where j corresponds to 
each of the three branches. In figl2](b) we present the result of adding pairs of branches, 
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N=1 




FIG. 2: (a) N = 1 BOMCA approximation corresponding to each of the three branches j = 1,2,3 
(see figCE]), \ipj(x,tf)\ = | exp[iSj(x, tf)/h]\. In the region of the transmitted wavefunction (starting 
from x ~ —0.2), a good approximation to the exact result is obtained by just a single branch 
(branch (1)). (b) The result of adding contributions from pairs of branches = [fy + t/ji\, i ^ j- 
Different choices of i and j are useful in different regions of the wavefunction. 

IV'I = l^j + V'il) * 7^ J- We see that to a good approximation the rippled parts of the 
wavefunction are described by a simple superposition with a proper choice of i and j. 

Note that as we approach the region of the transmitted part of the wavefunction (x ^ 
—0.2) a single branch (branch (1)) is enough to provide a good approximation. In fact, this 
is the branch that is responsible for the transmitted wavefunction (xf > 0, not shown), and 
hence we will call this the "transmitted branch". At energies on the order of magnitude of 
the barrier height, the "transmitted branch" is the real branch, but at lower energies (and 
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longer time scales) there is a crossover and the transmitted branch is one of the secondary 
branches. The existence of a single transmitted branch coincides well with the nodeless 
character of the transmitted part of the wavefnnction. By the sanre tohen, in Re, Q we 
showed that the reflected part of the wavefunction can be well approximated by a single 
branch if it has no ripples or nodes. 

The issue of which branches should be included in the sum and when, is obviously of 
central importance to the method; at present, we do not have a rigorous justification for 



the neglect of certain branches. A partial discussion is given in Ref. [21( in the context of 
GGWPD, but the BOMCA formulation, being more general, requires a more comprehensive 
discussion, which we leave to a future publication. 

We now turn to the next order of approximation, N = 2. The equations of motion are 
based on eqs. ffTUl) with two changes. The second equation in ffTQ]) is replaced with eq.©, in 
which there is a quantum force. Knowledge of the quantum force requires an equation of 
motion for v xx . We derive such an equation by inserting n = 2 in eqs. (ED and taking v {i) = 0, 
leading to 

dv X x Vxxx o /1 1 \ 

—7- = 3v x v xx . (11) 

at m 

At the level of iV = 2, four branches contribute to the reflected wavefunction at tf. Figj3]^a) 
depicts the contribution of the four branches, j = 1, ..,4. The result of adding pairs 
of branches = + i ^ j is given in figj3]^b). There is a significant increase in 
accuracy relative to iV = 1 in the vicinity of the maximum and to its right, although there 
is a slight decrease in accuracy relative to N = 1 to the left of the maximum. 

In conclusion, we have demonstrated that BOMCA accounts for interference and nodal 
structures of wavefunctions in a simple and natural way. In spite of the conceptual difficulties 
that crossing trajectories may pose for BM as an interpretational tool of QM, this notion 
introduces a powerful numerical tool and might even enrich the orthodox interpretation 
of the Bohmian formulation. The present results, combined with the simple and accurate 
results obtained for tunneling in Ref. [14J], demonstrates great promise for BOMCA as a 
versatile alternative to current semiclassical methods. As mentioned above, several issues 
require a more comprehensive understanding. First, what are the convergence properties of 
the method as higher order approximations are taken to the quantum force, that is increasing 
the value of iV? Can rigorous rules be derived for the summation of the different branches? 
What is the relation between the exact phase that diverges at a node and the approximate 
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N=2 




FIG. 3: (a) N = 2 BOMCA approximation corresponding to four branches \ipj(x,tf)\ = 
\exp[iSj(x,tf)/h]\, j = 1,...,4. (b) The result of adding contributions from pairs of branches 
I "01 = IV'j + ^i\-> i 7^ j- Except for a slight decrease in accuracy left to the maximum, the results 
for N = 2 where a quantum force term is added are better than for N = 1 (complex classical 
trajectories). 

BOMCA formulations that can account for nodes via a bipolar or multipolar expansion? 
We intend to address all these questions in future work. This work was supported by the 
Israel Science Foundation (576/04). 
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